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An Ising-like model of proteins is used to investigate the mechanical unfolding of the Green Flu- 
orescent Protein along different directions. When the protein is pulled from its ends, we recover the 
major and minor unfolding pathways observed in experiments. Upon varying the pulling direction, 
we find the correct order of magnitude and ranking of the unfolding forces. Exploiting the direction 
dependence of the unfolding force at equilibrium, we propose a force sensor whose luminescence 
depends on the applied force. 

PACS numbers: 87.15.A-; 87.15.Cc; 87.15.La 



INTRODUCTION 

In recent years many efforts have been devoted to the 
study of the mechanical properties of biopolymers un- 
der mechanical loading. Many experimental groups have 
studied the unfolding and refolding trajectories of pro- 
teins and nucleic acids by applying a controlled force 
with AFM or optical tweezers techniques [iHa]. These 
works have triggered a number of numerical investiga- 
tions, where the same molecules have been studied, under 
conditions otherwise not accessible to the experimental 
techniques [6l4l7 . 

One of the most interesting proteins studied with force 
spectroscopy is the Green Fluorescent protein (GPP), 
which exhibits bright green fluorescence when exposed 
to light with a suitable wavelength. GPP has many ap- 
plications in biotechnology, from localization of proteins 
in a living cell, to metal ion or pH sensor [18[. Purther- 
more, such a molecule has been the subject of mechanical 
experiments and numerical simulations [3|, |j, |8| aimed 
at characterizing its response to external force and the 
structure of its intermediate states. The final goal of such 
studies is a full characterization of the GPP response to 
mechanical stress, so as to pave the way to its use as a 
molecular force sensor. Indeed, it is commonly believed 
that GPP fluoresces only when its structure is almost 



(non-fluorescent) protein indicates that the applied force 
is below (above) some typical rupture force, but one can- 
not obtain an estimate of the actual value of the force by 
measuring the fluorescence. Thus in the present letter we 
propose a practical method to circumvent this limitation, 
by exploiting the fact that if pulled along different direc- 
tions, the GPP exhibits different mechanical properties, 
and thus different rupture forces, as already observed in 
experiments [4|. We will first introduce a model for the 
GPP that has already been used to evaluate the phase 
diagram, the free energy landscape |9| and the unfold- 



ing pathways [13 
RNA molecules 



mtact 



This represents a restriction for the use 
of the GPP to probe forces in vivo, since a fiuorescent 



J16| of widely studied proteins and of 
15|. We will compare the response of 
the model protein to experimental outcomes, and then 
we will study the rupture force of the GPP when pulled 
along different directions. Finally, we will introduce a 
model polyprotein made up of different GPP modules, 
and we will show how such a molecule can easily provide 
the value of the applied force in a wide range of values, 
and thus be used as a force probe. To the best of our 
knowledge, it is the first time such a molecular device is 
proposed in the literature. 



THE MODEL 

A native-centric, Ising-like model of protein folding 
23| has been generalized in previous works [9[ to deal 
with the case of mechanical unfolding. In such a model 



the state of a A^ residues long protein is determined 
by two sets of binary variables: the variables nik are 
associated to each residue being 1 or according to 
whether the residue is native-like or not, and the vari- 
ables aij = ±1 give the orientation (parallel or antiparal- 
lel to the external force) of a native-like stretch delimited 
by residues i and j > z in a non-native state, such that 
Sij = ii—'mi){l—mj)Y['l^i^i'rnk- A negative interaction 
energy hij is associated to two residues i and j if they 
are in contact in the native structure and if they are in 
the same native stretch. The corresponding hamiltonian 
is 

Af-l N j 

H{{m},{a})^Yl E h.jU'^k + UiLilm},^)), 



=1 j=i+i 



(1) 



b-l 



hjCTijSij is the molecule elongation 



where L — y > 

i—a j—i-\-l 

(distance between the Cq atoms of the two residues a 
and b where force is applied, 1 < a < b < N), lij is the 
native distance between Cq atoms of residues i and j , and 
C/(L({m} , {cr})) is the term describing the coupling to 
the external force, which depends on the pulling protocol. 
In the case of a constant force the energy reads U = —fL. 
In the case of a constant velocity protocol the force is 
applied through a harmonic potential whose center moves 
at constant velocity, and the corresponding energy is [/ = 
k/2{L-vt)^. 

Under a constant force the equilibrium thermodynam- 
ics of the model is exactly solvable [9[ , while for the study 
of the kinetics we resort to Monte Carlo simulations. 

Before studying the molecule response to external 
forces, it is interesting to consider the free energy pro- 
file at zero force. More specifically, we have computed 
the free energy profile as a function of the fraction of na- 
tive residues M (Fig. [Ij and contacts Q (Fig. [2]) at the 
denaturation temperature. Inspection of these plots in- 
dicates that at this temperature: (i) when the protein is 
in its native state, all the native contacts are formed, and 
almost all the residues are in the native configuration; (ii) 
in the unfolded state, no native contacts are formed, and 
1/3 of the residues are in the native configuration; (iii) 
the transition state corresponds to Q ~ 0. 3-^0. 4, while 
at Q -^ 0.5 ^0.7 there are some hints of the possible exis- 
tence of intermediate states; (iv) the unfolding barrier is 
of the order of 25 ksT in both cases. Our results for the 
free energy profile as a function of Q can be compared to 
the result obtained in [ll| by weighted-histogram anal- 
ysis of molecular dynamics (MD) data. The qualitative 
picture is very similar, although some differences can be 
observed. The MD result show that some fluctuations in 
native contacts are allowed in both the native and un- 
folded states, a feature which is missing in our result due 
to the extreme coopcrativity of the model. Moreover, 
the unfolding barrier is predicted by MD to be around 



15 fcsT: this is consistent with the observation that our 
model predicts systematically higher energy barriers and 
unfolding forces (see discussion below) . 







Figure 1. Free energy profile AG, as a function of the fraction 
of native residues M, at unfolding temperature T = 356 K 
and zero force. 
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Figure 2. Free energy profile AG, as a function of the fraction 
of native contacts Q, at unfolding temperature T = 356 K and 
zero force. 

From now on we set the temperature at T = 293 K. 
The native structure of GFP is basically a /3-barrel made 
of 11 /3-strands, with a N-terminal a-helix. 

In Fig. [3] the free energy profile AG{L) is reported 
for a typical case, where the equilibrium unfolding force 
/ = 35.9 pN is applied to the molecule ends. Besides the 
native and the unfolded minima we can see three other 
local minima (or bends which become actual minima at 
different values of the force) around 11, 18 and 25 nm. 
As will be shown in detail in the next section, these lo- 
cal minima and bends correspond to intermediate states 
effectively populated in the simulations. Analysing the 
equilibrium probability < (mfc(L)) < 1 that the fc— th 
residue is native-like when the molecule total elongation 
is L (data not shown), wc find out that such bends corre- 
spond to the following structures: f3i and /3ii (for i ~ 11 



nm), /3io/3ii (18 nm) and /3i;52/33 (25 nm). Here and 
in the following, /3k ■ ■ ■ Pn denotes an unfolded structure 
of the GFP, where /3~strands from fc to n are not in a 
native-like conformation, i.e. they are unfolded (in all 
these structures the N-terminal a-helix is also not in a 
native-like conformation) . 







Figure 3. Free energy profile AG as a function of the protein 
elongation L, with T = 293 K, and force / = 35.9 pN applied 
to the molecule ends. 



SIMULATIONS AND COMPARISON WITH 
EXPERIMENTS 

We first consider simulations at constant velocity, 
mimicking the effect of an AFM cantilever, which is re- 
tracted at a velocity v. The force is applied to the 
molecule ends and we set /c = 30 pN/nm and consider 
velocities v = 0.3, 1 /xm/s, 2 /xm/s and 3.6 /xm/s. In 
ref. [8| it was found that the most likely unfolding path- 
way is a — > /?i — >^ /52/33 (the a-hclix unfolds first, then 
/3-strand 1, then strands 2 and 3 simultaneously, then 
the rest of the protein), observed in 72% of the trajecto- 
ries at w = 2.5 ^m/s. A different pathway, a — >■ /3ii, was 
observed in the remaining 28% of the trajectories. 

The N-terminal a-helix is the first secondary struc- 
ture element to unravel in both pathways. This event is 
typically associated to very small signals, almost masked 
by fluctuations, at odds with the clear jumps we observe 
in the end-to-end length for the detaching of /3-strands 
(see fig. 2]). This is analogous to what occurs in the ex- 
periments where the unfolding of the helix is associated 
to a very smooth "hump-like" transition with a short 
contour length increase of 3.2 nm in the force-extension 
traces in [3| and by a small jump in the root mean square 
distance as a function of time in [8| . 

We have found that, at all velocities considered, in less 
than 10% of the trajectories /3ii is the first strand to 
unravel, while the remaining trajectories follow the ma- 
jor unfolding pathway found in experiments. In Fig. 21 



we plot a typical unfolding trajectory of our GFP model 
when the force is applied to the molecule ends. In partic- 
ular, we plot, as functions of time, the end-to-end length 
L, the force and several weighted fractions of contacts 
between adjacent /3-strands ^^^.^ , giving the fraction of 
native contacts between the strands /3j and Pj , see [l6| . 

Major unfolding pathway. Inspection of the top panel 
of Fig. m corresponding to the major unfolding pathway, 
provides clear evidence that there are three main unfold- 
ing events, (i) A drop in the number of contacts between 
strands 1 and 2, signalling the unfolding of /3i (actually 
the a-helix has already unfolded, as discussed above, but 
the corresponding fraction of native contacts 4>a is not re- 
ported in the figure for the sake of clarity). The length 
of the corresponding intermediate state is in the range 
10-^ 12.5 nm, where the free energy profile AG{L) shows 
a bend, see fig. [H (ii) A drop in the number of con- 
tacts involving strands 2 and 3, signalling the unfolding 
of these strands. The corresponding intermediate length 
is around 20 nm, where the free energy profile AG(L) has 
a local minimum, (iii) A drop in the number of contacts 
involving strands 10 and 11, signalling the unfolding of 
these strands. The corresponding intermediate length is 
in the range 30-f37 nm: inspection of fig. |3] suggests that 
for such an elongation the molecules lies in the basin of 
the unfolded minimum. 

Minor unfolding pathway. The bottom panel of Fig.|4] 
corresponds to the minor unfolding pathway and we can 
see that the first strand to unravel is /3ii followed by /3io. 
In [8| the unfolding pathway was only traced up to the 
A/3ii intermediate because the subsequent event is the 
fiattening of the barrel but, after the barrel flattens, there 
is at least another rupture event as the last force jump in 
Fig. lb of [8| shows. It is reasonable to assume that this 
event is related to the breaking of native-like contacts 
between the beta strands, which were not ruptured dur- 
ing the flattening of the barrel. Our model, which lacks 
a fully three-dimensional representation, cannot describe 
the flattening of the barrel, while it can describe with 
a high time resolution the breaking of the beta strand 
contacts, which here yield in a few distinct steps (Fig. |4j 
bottom panel). 

We can now put the local minima and bends of the 
free energy landscape of Fig. |3] (which is a thermody- 
namic equilibrium property of the system) in correspon- 
dence with intermediates found in our simulations and 
in experiments (which are performed in non-equilibrium 
conditions). Some of these features of the free energy pro- 
file are indeed barely visible, but the equilibrium proba- 
bilities (TOfc(L)), we introduced in the previous section to 
give a structural interpretation of the various minima and 
bends, are perfectly consistent with the nonequilibrium 
TOfc values obtained from the simulations, which allow us 
to identify the structures of the nonequilibrium interme- 
diates. 

In ref. [3[ the authors observed two intermediates with 
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Figure 4. (color online) Typical unfolding trajectories of a 
GFP module under constant velocity pulling {v = 0.3 /im/s). 
Length L, force / and weighted fractions (l>^i- of strand- 
strand contacts as functions of time for two typical cases: 
major (top panel) and minor (bottom panel) unfolding path- 
way, see text. 



separation values from native configuration of 3.2 and 
10 nm. The first one, is an intermediate with only the 
N-terminal a-helix detaehed that profile of Fig. [3] does 
not show, while the seeond is an intermediate with the 
N-terminal a-hefix detached and a /3-strand detached 
which corresponds to the bend at 11 nm (9.2 nm awa' 
from native state) in Fig. [S] The authors of ref. 
report the existence of another intermediate (N-terminal 
a~hcfix and first, second and third /3-strands detached) 
with a distance of 26.3 nm from the native statc(16.3 nm 
from the previous second intermediate) which is clearly 
associated to our dip at 25 nm. The 18 nm intermediate 
has no analogue in experiments. 

Different directions. We now consider simulations 
where the points of force application are not the molecule 
ends, so that the direction of the force with respect to 
the molecule is varied. Table U reports, for different di- 
rections (specified through the application point residue 
numbers), the mean unfolding forces, where unfolding is 
defined as unravelling of the first /3 strand. Since most 
of these directions were considered in experiments [J] , at 



least at u = 3.6 /im/s, it is interesting to compare our 
results to the experimental ones. Our unfolding forces 
are systematically larger than the experimental values, 
with the largest discrepancies (a factor 2 to 3) occur- 
ring for directions 3-212 and 132-212. However, it is 
interesting that in spite of the simplicity of the model, 
which lacks a fully three-dimensional representation, the 
orders of magnitude for the rupture forces are correct 
and many qualitative aspects are reproduced. In partic- 
ular, by analyzing the experimental data one finds that 
the unfolding force increases with the following order: (i) 
pulling along the end-to-end direction (it must be noted 
that the rupture force along this direction was measured 
for V = 0.3 ^m/s instead of 3.6 //m/s as most other di- 
rections); (ii) pulling along the 3-212 and 132-212 direc- 
tions, the corresponding rupture forces are equal within 
the experimental error; (iii) pulling along the 182-212 
and 3-132 directions, the corresponding rupture forces 
are equal within the experimental error (though the lat- 
ter was measured for v = 2 /xm/s) ; (iv) pulling along the 
117-182 direction. This hierarchy is respected by our re- 
sults: we find that the rupture force increases when we 
consider the pulling directions as ordered above, the only 
exception being for 3-212 and 132-212, whose unfolding 
forces are not equal (we obtain a smaller force for the 
latter), and the same holds for 182-212 and 3-132 (we 
obtain a larger force for the latter). 



Table I. Unfolding forces at different velocities for different 
directions. Experimental values (* from ref. [3] and ^ from 
ref. [J]) in parentheses. 



Direction 


V = 0.3 /im/f 


Unfolding force 
i V = 2 /im/s 


(PN) 

V = 3.6 /im/s 


end-end 


140 ±3 
(104 ±40)* 


177 ±7 


184 ± 13 


182-end 


196 ±7 


226 ±6 


244 ±7 


3-212 


244 ± 12 


298 ± 12 


317 ±20 
(117±19)t 


132-212 


251 ±7 


266 ±3 


273 ±6 

(127±23)t 


132-end 


306 ± 12 


360 ± 20 


381 ± 26 


182-212 


365 ±2 


390 ±7 


409 ± 15 
(356±61)t 


3-132 


383 ± 16 


471 ± 49 
(346 ± 46) t 


535 ± 80 


117-182 


467 ±3 


501 ±11 


512 ±11 

(548±57)t 



Finally, in TablelUwe report the potential width values 
Xu corresponding to the rupture of the first /3-strand for 
different directions. These were obtained through a fit 
of the most probable unfolding force /m as a function of 



velocity to the Evans-Ritchie theory [2l|, [22] , which gives 



/ 



M 






(2) 



where tq is the unfolding time at zero force. It must be 
kept in mind that in the Evans-Ritchie theory the force 
grows with a constant rate r = k ■ v and hence its ap- 
plicability to the present case (harmonic potential whose 
center moves at constant velocity v) is only approximate. 



Table II. Potential width Xu obtained from a fit to Eq. (2] 
Experimental values between parentheses. 



direction 


Xu (nm) 


direction 


Xu (nm) 


end-end 


0.21 ±0.04 
(0.28 ± 0.03) 


132-end 


0.14 ±0.01 


182-end 


0.22 ±0.03 


182-212 


0.24 ±0.04 
(0.14 ±0.002) 


3-212 


0.14 ±0.01 
(0.45 ±0.01) 


3-132 


0.11 ±0.03 
(0.125 ±0.005) 


132-212 


0.46 ± 0.06 
(0.32 ± 0.005) 


117-182 


0.22 ±0.01 
(0.12 ±0.003) 




Figure 5. (color online) Sketch of a polyprotein made of 
various modules connected between them through different 
residues. 



made up of several GFP modules, where a module linked 
by its (3,212) residues to the main structure was alter- 
nated with a module linked by its (132,212) residues. 
Such a molecule can be easily obtained by using the cys- 
teine engineering method discussed in ref. [26| , which al- 
lows one to construct polyproteins with precisely con- 
trolled linkage topologies: the points of force application 
to each module correspond to the position of the linking 
cysteines in the folded tertiary structure. In order to un- 
derstand the general behavior of our model polyprotein 
under a constant force we first investigate the response 
to a constant force of a single GFP module. The corre- 
sponding unfolding forces are reported in Table [ 



Our potential widths are consistent with experimental 
ones only in a few cases (end-end, 3-132): once again, 
this might be attributed to the fact that our model lacks 
a fully three-dimensional representation. 

It must also be observed that the Evans-Ritchie theory 
is built on the assumption that Xu is independent of the 
applied force, and this can be another source of error in 
the determination of x,,, ■ This assumption was relaxed in 
more recent theories [2J, |25J which yield generalizations 
of Eq. [21 which predict that the /m vs In v plot is non- 
linear, with the slope being an increasing function of v, 
as observed in many experiments. Indeed our data show 
some nonlinearity, but this is too small to apply these 
theories, probably because our velocities span only 1 or- 
der of magnitude. Previous applications of these theories 
[9|,|2J,|25| were done on data sets with velocities spanning 
4-5 orders of magnitude, such that the nonlinear effects 
were much more important, but such a wide range of 
velocities is beyond the scope of the present paper. 



GFP POLYPROTEIN 

As we discussed above, the equilibrium properties of 
the GFP at constant force, can be obtained exactly, for 
any pulling direction. Exploiting this result, we study a 
polyprotein where each module is connected to the neigh- 
bouring ones through different points of force application, 
as illustrated in Fig. [S] 

For example Dietz et al. [4| already proposed a copoly- 
mer with mixed linkage geometries GFP(3, 212) (132,212), 



Table III. 


Equilibrium unfolding force for different directions. 


direction unfolding force 


(pN) 


direction unfolding force (pN) 


end-end 




35.9 




182-end 


65.0 


3-212 




38.9 




117-190 


67.3 


3-182 




42.6 




102-190 


71.2 


117-end 




50.8 




117-182 


78.1 


3-132 




56.4 




132-182 


96.7 



We then proceed by studying the response to a con- 
stant force of a polyprotein made up of 10 modules, each 
with different linkage topologies. It is worth to note that 
at equilibrium a force applied to the free ends of the 
polyproteins will have the same value throughout the 
whole chain. Thus, the different modules will unfold at 
different values of the force, according to the hierarchy 
shown in Tab. IIIIl and thus the luminescence will be dif- 
ferent for different values of the force. If we assign a 
value 1 (in an arbitrary scale) to the maximum possible 
luminescence, where each module is emitting green light, 
a luminescence of 0.5 will correspond to a configuration, 
and thus to a force, where half of the modules are un- 
folded (non intact structure). Given that each module 
with a different linkage has a different unfolding force, 
we obtain a curve like the one shown in fig. [6l relating 
the luminescence of the polyprotein to the force applied 
to its free ends, where the force ranges from 35.9 to 96.7 
pN, see Table IIIIl It is worth to note that interface in- 
teractions and aggregation effects between neighbouring 
units in polyproteins similar to the one we propose, have 



been ruled out by experimental investigations, see 
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Figure 6. Fraction of native-like modules as a function of 
force at T =293 K. Each "step" corresponds to the unfolding 
of a different module in the polyprotein and thus to a decrease 
in the luminescence by a "unit" . 



dimensional representation. Moreover, from a more 
quantitative point of view, our energy barriers and un- 
folding forces arc systematically larger than those ob- 
served in experiments. 

We have exploit the dependence of the unfolding force 
on the pulling direction to investigate a force sensor based 
on a GFP polyprotein where each module is linked with 
a different geometry to the nearest neighbouring mod- 
ules, so as to experience the force along different direc- 
tion, yielding a device whose luminescence depends (in 
a discrete way) on the force. It is worth noting that 
such a device may be used in in-vivo experiments, to 
measure forces at molecular level, e.g. inside cells, in a 
non-invasive way. 

MC gratefully acknowledges the financial support of 
Aarhus Universitets Forskningsfond (AUFF) during his 
stay at Aarhus University. AI gratefully acknowledges fi- 
nancial support from Lundbeck Fonden. AP thanks An- 
tonio Trovato for a useful discussion. 



Figure [6] represents an important result of this paper. 
It is worth to note that more modules, with different link- 
ages, can be added, and this would give a more precise 
determination of the force. Once the polyprotein we pro- 
pose has been engineered, a curve like the one shown in 
fig. |6] can be very easily obtained in an optical tweezers 
experiment at constant force as those discussed, e.g., in 
ref. PJ. This approach would also allow one to calibrate 
the device. Finally, we want to emphasize that, although 
unfolding studies of GFP alongdifferent directions where 
already performed in, e.g., |j, ISfj those previous studies 
considered the dynamic-loading set up, with a constant 
retraction speed of the AFM cantilever. On the contrary 
we investigate here for the first time the unfolding at 
constant force of GFP. The unfolding force of a molecule 
under dynamic loading depends not only on the molecu- 
lar features, but also on the force rate, and thus a force 
probe based on those data must be able to measure at the 
same time the loading rate and the rupture force. Our 
constant force probe does not exhibit this drawback. 



CONCLUSIONS 

The Ising-like model of protein mechanical unfolding 
describes correctly the most important qualitative as- 
pects of the direction-dependent mechanical unfolding 
of the Green Fluorescent Protein, namely the unfold- 
ing pathways and intermediates observed when pulling 
at constant velocity from the molecule ends, and the or- 
ders of magnitude and ranking of the unfolding forces 
corresponding to different directions. Some features, 
like the flattening of the barrel or the potential widths 
corresponding to many directions, cannot however be 
described by our model, which lacks a fully three- 
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